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Abstract 

A residual-based  lubrication  method  is  used  in  this  paper  to  find  the  flow  rate 
and  pressure  field  in  converging-diverging  rigid  tubes  for  the  flow  of  time- 
independent  category  of  non-Newtonian  fluids.  Five  converging-diverging 
prototype  geometries  were  used  in  this  investigation  in  conjunction  with 
two  fluid  models:  Ellis  and  Herschel-Bulkley.  The  method  was  validated  by 
convergence  behavior  sensibility  tests,  convergence  to  analytical  solutions  for 
the  straight  tubes  as  special  cases  for  the  converging-diverging  tubes,  conver- 
gence to  analytical  solutions  found  earlier  for  the  flow  in  converging-diverging 
tubes  of  Newtonian  fluids  as  special  cases  for  non-Newtonian,  and  conver- 
gence to  analytical  solutions  found  earlier  for  the  flow  of  power-law  fluids 
in  converging-diverging  tubes.  A brief  investigation  was  also  conducted  on 
a sample  of  diverging-converging  geometries.  The  method  can  in  principle 
be  extended  to  the  flow  of  viscoelastic  and  thixotropic/rheopectic  fluid  cat- 
egories. The  method  can  also  be  extended  to  geometries  varying  in  size  and 
shape  in  the  flow  direction,  other  than  the  perfect  cylindrically-symmetric 
converging-diverging  ones,  as  long  as  characteristic  flow  relations  correlat- 
ing the  flow  rate  to  the  pressure  drop  on  the  discretized  elements  of  the 
lubrication  approximation  can  be  found.  These  relations  can  be  analytical, 
empirical  and  even  numerical  and  hence  the  method  has  a wide  applicability 
range. 

Keywords:  non-Newtonian  fluids;  Ellis;  Herschel-Bulkley;  yield-stress;  power- 
law;  Bingham;  converging-diverging  tubes;  diverging-converging  tubes;  non- 
linear systems. 


1 Introduction 

The  flow  of  fluids  in  general  and  non-Newtonian  in  particular  through  conduits 
with  varying  cross  sectional  size  and/or  shape  is  commonplace  in  natural  and  in- 
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dustrial  flow  systems  such  as  filtration  and  refinement  devices  in  the  chemical  in- 
dustries and  biological  fluid  transportation  networks  like  blood  circulation  vessels 
and  air  respiration  pathways  [1-7].  Modeling  such  flows  is  also  needed  for  analyz- 
ing fluid  transportation  through  porous  media  where  the  irregularly-shaped  pores 
and  throats  are  usually  described  by  geometrically-simple  flow  conduits  with  vary- 
ing cross  sections  in  the  flow  direction  [8-17].  In  this  respect,  pore-scale  network 
modeling  is  the  main  candidate  and  the  most  natural  approach  to  accommodate 
such  flow-structure  features  [18-23].  The  effect  of  varying  cross  section  in  the  flow 
direction  can  be  particularly  important  for  modeling  certain  flow  phenomena  such 
as  extensional  flows,  viscoelastic  effects  and  yield-stress  dynamics  which  are  essen- 
tial for  various  scientific  and  technological  purposes  like  enhanced  oil  recovery  as 
well  as  many  other  applications  [24-35]. 

Several  methods  have  been  used  in  the  past  to  model  and  simulate  such  flows; 
these  include  numerical  methods  like  finite  element,  finite  difference  and  stochastic 
algorithms  [9-11,  36-38],  as  well  as  analytically-based  methods  [1,  2,  39,  40].  Most 
previous  attempts  employ  complex  mathematical  and  computational  techniques, 
which  may  be  unwanted  due  to  difficulties  in  implementation  and  verification,  high 
computational  costs,  computational  complexities  like  convergence  difficulties,  as 
well  as  susceptibility  to  errors  and  bugs.  Also,  many  of  the  previous  investigations 
are  relevant  only  to  special  cases  of  fluid  or  flow  or  conduit  geometry  with  a sub- 
sequent limited  use.  Therefore,  developing  a simple,  general  and  robust  method  to 
investigate  such  flows  is  highly  desired. 

In  this  paper  we  use  a residual-based  lubrication  approximation  method  to  find 
the  flow  of  non-Newtonian  fluids  through  rigid  tubes  with  converging-diverging 
and  diverging-converging  shapes  in  the  axial  direction.  The  method  is  based  on 
discretizing  the  flow  conduits  in  the  axial  dimension  into  ring-like  elements  on 
which  non-Newtonian  characteristic  flow  relations  derived  for  conduits  with  a fixed 
axial  geometry  can  apply.  The  formulation  in  this  paper  is  limited  to  the  time- 
independent  category  of  the  non-Newtonian  fluids;  although  the  method  in  prin- 
ciple is  capable  of  being  extended  to  history-dependent  fluids,  i.e.  viscoelastic 
[24]  and  thixotropic/rheopectic.  Two  widely  used  non-Newtonian  fluid  models, 
Ellis  and  Herschel-Bulkley,  are  examined  in  this  study  as  prototypes  for  the  time- 
independent  models  to  which  this  method  applies.  These  two  fluids  can  be  used  to 
model  diverse  non-Newtonian  rheological  phenomena  such  as  shear-thinning,  shear- 
thickening, and  yield  stress  as  well  as  Newtonian  flow  as  a special  case.  These  fluid 
models  are  commonly  used  for  describing  the  time-independent  non-Newtonian 
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rheology  especially  in  the  industrial  applications  such  as  polymer  manufacturing 
and  processing  and  oil  recovery  [41,  42], 

Although  the  present  paper  investigates  the  flow  of  non-Newtonian  fluids  in  cer- 
tain axi-symmetric  geometries  with  converging  and  diverging  features,  the  method 
is  more  general  and  can  be  used  with  other  geometries  whose  cross  section  varies 
in  size  and/or  shape  in  the  flow  direction.  The  only  condition  for  the  applicability 
of  the  method  is  the  availability  of  analytic,  empirical  or  numerical  [43]  relations 
that  correlate  the  flow  rate  in  the  discretized  elements  to  the  pressure  drop  across 
these  elements. 


2 Investigated  non-Newtonian  Fluids 

In  this  section  we  present  a general  background  about  the  bulk  rheology  and  flow 
relations  for  the  two  investigated  non-Newtonian  fluids:  Ellis  and  Herschcl-Bulkley. 

2.1  Ellis  Model 

This  is  a three-parameter  model  that  is  widely  employed  to  describe  the  flow  of 
time-independent  shear-thinning  yield-free  non-Newtonian  fluids.  The  model  is 
known  for  being  successful  in  describing  the  bulk  and  in  situ  rheologies  of  a wide 
range  of  polymeric  solutions.  It  is  particularly  used  as  a replacement  for  the  power- 
law  model  due  to  its  superiority  in  matching  experimental  data.  The  model  is 
characterized  by  a single  low-shear-rate  Newtonian  plateau  with  the  absence  of  a 
high-shear-rate  one  and  hence  it  may  better  match  the  observations  in  low  and 
medium  shear-rate  regimes  than  in  the  high  shear-rate  regimes. 

An  attractive  feature  of  the  Ellis  model  is  the  availability  of  an  analytical  ex- 
pression that  links  the  flow  rate  to  the  pressure  drop  for  the  flow  in  cylindrical 
constant-radius  tubes.  This  feature  makes  rheological  modeling  with  Ellis  more 
accurate  and  convenient  than  with  some  other  models,  such  as  Carreau,  which 
have  no  such  analytical  expressions  although  these  models  may  be  as  good  as 
or  even  superior  to  Ellis.  The  advantage  of  having  a closed-form  analytical  flow 
characterization  relation  is  obvious  especially  for  some  numerical  techniques  like 
pore-scale  network  modeling. 

According  to  the  Ellis  model,  the  fluid  viscosity  p as  a function  of  the  shear- 
stress  is  given  by  the  following  expression  [21,  44-48] 


3 


a— l 


(i) 


11 


t^G 


1 + I 7^ 

Tl/2 


where  /iQ  is  the  zero-shear-rate  viscosity,  r is  the  shear  stress,  r1/2  is  the  shear  stress 
at  which  p,  — and  ct  is  a flow  behavior  index  in  this  model. 

For  Ellis  fluids,  the  volumetric  flow  rate  as  a function  of  the  pressure  drop  in  a 
circular  cylindrical  tube  with  a constant  radius  over  its  length  assuming  a laminar 
incompressible  flow  with  no  slip  at  the  tube  wall  [49]  is  given  by  the  following 
equation  [21,  50] 


nR^Ap  4 / RApY  1 

8 Lp0  +a  + 3\^2Lr1/2y 


(2) 


where  R is  the  tube  radius,  A p is  the  pressure  drop  over  the  tube  length  and  L is 
the  tube  length.  The  derivation  of  this  expression  is  given  in  [21,  50]. 


2.2  Herschel-Bulkley  Model 

This  is  a three-parameter  model  that  can  describe  Newtonian  and  large  group  of 
time-independent  non-Newtonian  fluids.  According  to  the  Herschel-Bulkley  model, 
the  shear  stress  as  a function  of  the  shear  rate  is  given  by  the  following  relation 
[21,  48,  51] 

r = t0  + C'y71  (r  > r0)  (3) 

where  r is  the  shear  stress,  tq  is  the  yield-stress  above  which  the  substance  starts  to 
flow,  C is  the  consistency  coefficient,  7 is  the  shear  rate  and  n is  the  flow  behavior 
index.  The  Herschel-Bulkley  model  reduces  to  the  power-law,  or  Ostwald-de  Waele 
model,  when  the  yield-stress  is  zero,  and  to  the  Bingham  plastic  model  when  the 
flow  behavior  index  is  unity.  It  also  emulates  the  Newtonian  fluids  when  both  the 
yield-stress  is  zero  and  the  flow  behavior  index  is  unity  [52]. 

The  Herschel-Bulkley  model  contains  six  classes:  (a)  shear-thinning  without 
yield-stress  [n  < 1.0,  r0  = 0]  which  is  the  power-law  fluid  in  shear-thinning  mode, 
(b)  shear-thinning  with  yield-stress  [n  < 1.0,  r0  > 0],  (c)  Newtonian  [n  = 1.0,  ra  = 
0],  (d)  Bingham  plastic  [n  = 1.0,  rD  > 0],  (e)  shear-thickening  (dilatant)  with- 
out yield-stress  [n  > 1.0,  rQ  = 0],  and  (f)  shear-thickening  with  yield-stress  [n  > 
1.0,  t0  > 0], 

For  Herschel-Bulkley  fluids,  the  volumetric  flow  rate  as  a function  of  the  pres- 
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sure  drop  in  a circular  cylindrical  tube  with  a fixed  radius  assuming  a laminar 
incompressible  flow  with  no  wall  slip  [49,  53]  is  given  by  the  following  relation 
[21,  50,  51]: 


Q 

Q 
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2 Tp  (tw  - T0) 

2 + 1 /n 
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where  tw  is  the  shear  stress  at  the  tube  wall  which  is  given  by  tw  = while 
the  meaning  of  the  other  symbols  have  already  been  given.  The  derivation  of  this 
expression  can  be  found  in  [21,  50]. 


3 Method 

According  to  the  residual-based  lubrication  approach,  which  is  proposed  in  the 
present  paper  to  find  the  flow  rate  and  pressure  field  in  converging-diverging  tubes 
for  the  flow  of  time- independent  non-Newtonian  fluids,  the  tube  is  discretized  in  the 
axial  dimension  into  ring-like  elements.  Each  one  of  these  elements  is  treated  as  a 
constant-radius  single  tube,  where  the  radius  of  the  element  is  taken  as  the  average 
of  its  inlet  and  outlet  radii,  to  which  Equations  2 and  4 apply.  From  this  discretized 
form  of  the  flow  conduit,  a system  of  non-linear  simultaneous  equations,  which  are 
based  on  the  mass  conservation  residual  plus  boundary  conditions,  is  formed. 

For  a discretized  tube  consisting  of  (IV  — 1)  elements,  there  are  N nodes:  (IV  — 2) 
internal  nodes  and  two  boundary  nodes.  Each  one  of  these  nodes  is  characterized 
by  a well-defined  axial  pressure  according  to  the  one-dimensional  flow  model.  Also 
for  the  internal  nodes,  and  because  of  the  assumption  of  flow  incompressibility, 
the  total  sum  of  the  volumetric  flow  rate,  which  is  signed  (+/— ) according  to  its 
direction  as  being  toward  or  away  from  the  given  node,  is  zero  due  to  the  lack  of 
sources  and  sinks  at  the  node,  and  hence  (N  — 2)  residual  functions  whose  essence 
is  the  vanishing  of  the  net  flow  at  the  internal  nodes  are  formed.  These  residual 
equations  are  coupled  with  two  given  boundary  conditions  at  the  inlet  and  outlet 
nodes  to  form  a system  of  N simultaneous  equations. 

A typical  method  for  solving  such  a system  is  to  use  an  iterative  non-linear 
solution  scheme  such  as  Newton-Raphson  method  where  an  initial  guess  for  the 
internal  nodal  pressures  is  proposed  and  used  in  conjunction  with  the  Jacobian 
matrix  of  the  residual  functions  to  find  the  pressure  perturbation  vector  which  is 
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then  used  to  adjust  the  values  of  the  internal  nodal  pressure.  This  iterative  process 
is  repeated  until  a convergence  condition,  which  is  usually  based  on  the  size  of  the 
residual  norm,  is  reached. 

In  more  formal  terms,  the  process  is  based  on  solving  the  following  matrix 
equation  iteratively 


JAp  = — r (5) 

where  J is  the  Jacobian  matrix,  p is  the  vector  of  variables  which  represent  the 
pressure  values  at  the  boundary  and  internal  nodes,  and  r is  the  vector  of  residuals 
which,  for  the  internal  nodes,  is  based  on  the  continuity  of  the  volumetric  flow  rate 
as  given  by 


fj  = Y.&  = ° (6) 

i=l 

where  m is  the  number  of  discretized  elements  connected  to  node  j which  is  two 
in  this  case,  and  Qi  is  the  signed  volumetric  flow  rate  in  element  i as  characterized 
by  Equations  2 and  4.  Equation  5 is  then  solved  in  each  iteration  for  Ap  which 
is  then  used  to  update  p.  The  convergence  will  be  announced  when  the  norm  of 
the  residual  vector,  r,  becomes  within  a predefined  tolerated  marginal  error.  More 
details  about  this  solution  scheme  can  be  found  in  [37,  54,  55]. 
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4 Implementation  and  Results 


The  residual-based  lubrication  method  was  implemented  in  a computer  code  and 
flow  simulation  results  were  obtained  for  a wide  range  of  fluid,  flow  and  tube 
parameters  for  both  Ellis  and  Herschel-Bulkley  fluids.  All  the  simulations  reported 
in  the  present  paper  were  carried  out  using  evenly-divided  discretization  meshes. 
In  this  investigation  we  used  five  cylindrically-symmetric  converging-diverging  tube 
geometries,  a graphic  demonstration  of  which  is  shown  in  Figure  1.  The  equations 
that  describe  the  radius  dependency  on  the  tube  axial  coordinate  for  these  five 
geometries  are  given  in  Table  1,  while  a generic  converging-diverging  tube  profile, 
demonstrating  the  coordinate  system  used  in  the  mathematical  formulation  of  this 
dependency,  is  shown  in  Figure  2. 

A sample  of  the  results  of  these  simulations  are  shown  in  Figure  3 for  Ellis 
fluids  and  Figure  4 for  Herschel-Bulkley  fluids  using  the  five  geometries  of  Table 
1.  In  these  figures  the  tube  axial  pressure  is  plotted  as  a function  of  the  tube 
axial  coordinate  for  a number  of  authentic  fluids  whose  bulk  rheologies  are  given  in 
Table  2 for  the  Ellis  fluids  and  Table  3 for  the  Herschel-Bulkley  fluids.  The  pressure 
dependency  for  the  Poiscuille  flow  of  Newtonian  fluids  having  the  same  zero-shear- 
rate  viscosity  is  also  given  in  these  figures  for  demonstration  and  comparison.  The 
inclusion  of  the  Poiseuille  flow  can  also  serve  as  a sensibility  test;  for  example  in 
Figure  3 (c)  we  see  a strong  similarity  in  the  pressure  field  and  flow  rate  between  the 
Poiscuille  and  Ellis  flows  since  the  latter  has  a strong  resemblance  to  the  Newtonian 
at  these  flow  regimes  as  can  be  seen  from  its  rheological  parameters.  Also  for  the 
flow  represented  by  Figure  4 (e)  we  observe  the  convergence  of  the  Bingham  flow 
rate  to  the  Newtonian  flow  rate  at  this  high-pressure  flow  regime  which  is  a sensible 
trend.  Unlike  the  flow  rates  of  the  other  figures,  the  non-Newtonian  flow  in  Figure  4 
(e)  is  lower  than  the  Newtonian  flow  rate  because  the  non-Newtonian  is  a Bingham 
fluid  and  not  a shear-thinning  fluid.  The  effect  of  the  yield-stress  in  this  case  is  to 
lower  the  flow  rate  of  Bingham  from  its  Newtonian  counterpart,  as  expected.  It 
should  be  remarked  that  in  Figures  3 and  4,  pi  and  pQ  represent  the  inlet  and  outlet 
pressures,  while  Qe,  Qh  and  Qp  are  the  flow  rates  for  Ellis,  Herschel-Bulkley  and 
Poiseuille  respectively.  The  plots  shown  in  these  figures  represent  the  converged 
solutions  which  are  obtained  with  the  use  of  discretization  meshes  that  normally 
consist  of  50-100  elements. 

In  Figures  5-9  we  present  a sample  of  our  results  in  a different  form  where 
we  plot  the  volumetric  flow  rate  as  a function  of  the  pressure  drop  for  a num- 
ber of  converging-diverging  geometries  using  a representative  sample  of  Ellis  and 
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Herschel-Bulkley  fluids.  It  should  be  remarked  that  the  residual-based  lubrication 
method  can  also  be  used  with  diverging-converging  geometries  similar  to  the  ones 
demonstrated  in  Figure  10.  A sample  of  the  flow  simulation  results  using  the  lat- 
ter geometries  with  a representative  sample  of  Ellis  and  Herschel-Bulkley  fluids  is 
shown  in  Figures  11-13.  In  each  one  of  these  figures,  plots  of  both  axial  pressure 
versus  axial  coordinate  and  flow  rate  versus  pressure  drop  are  given.  An  interesting 
feature  is  that  the  curvature  of  the  pressure  field  plots  in  these  figures  is  qualita- 
tively different  from  that  seen  in  the  plots  of  the  converging-diverging  geometries, 
which  is  a sensible  tendency  as  it  reflects  the  nature  of  the  conduit  geometry. 


Table  1:  The  correlation  between  the  tube  radius  R and  its  axial  coordinate  x for 
the  five  converging-diverging  geometries  used  in  this  paper.  In  these  equations, 
< x < ^ and  Rm  < Rm  where  Rm  is  the  tube  minimum  radius  at  x = 0 and 
Rm  is  the  tube  maximum  radius  at  x = ±-|  as  demonstrated  in  Figure  2. 


Geometry 

Conic 

Parabolic 

Hyperbolic 


R(x) 

~D  | 2{RM—Rm)  I I 

1 1 in  I ^ | ^ | 

Rm  + (x)  (Rm  — Rm)%2 

\As,  + (if  (R2m  - Rl)*2 


Hyperbolic  Cosine 


Rm  cosh 


1-arccosh 


x 


Sinusoidal 


(RM+Rmj  _ cos  (2p) 


Table  2:  Bulk  rheology  data  related  to  the  Ellis  fluids  used  to  generate  the  plots 
in  Figure  3. 


Source 

Fluid 

Ho  (Pa.s) 

a 

Ti/ 2 (Pa) 

Fig.  3 (a) 

Sadowski  [56] 

0.4%  Natrosol  - 250H 

0.1000 

1.811 

2.2 

Fig.  3 (b) 

Sadowski  [56] 

1.4%  Natrosol  - 250G 

0.0688 

1.917 

59.9 

Fig.  3 (c) 

Sadowski  [56] 

6.0%  Elvanol  72-51 

0.1850 

2.400 

1025.0 

Fig.  3 (d) 

Park  [57] 

polyacrylamide  0.50% 

4.35213 

2.4712 

0.7185 

Fig.  3 (e) 

Park  [57] 

polyacrylamide  0.05% 

0.26026 

2.1902 

0.3390 
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(a)  Conic 


(b)  Parabolic 


(c)  Hyperbolic 


(e)  Sinusoidal 

Figure  1:  A graphic  demonstration  of  the  converging-diverging  tube  geometries 
used  in  this  investigation. 
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Figure  2:  A schematic  profile  of  the  converging-diverging  tubes  that  demonstrates 
the  setting  of  the  coordinate  system  used  in  the  equations  correlating  the  tube 
radius  R to  its  axial  coordinate  x as  used  in  Table  1. 


Table  3:  Bulk  rheology  data  related  to  the  Herschcl-Bulkley  fluids  used  to  generate 
the  plots  in  Figure  4. 


Source 

Fluid 

C (Pa.sn) 

n 

r0  (Pa) 

Fig.  4 (a) 

Park  [57] 

0.50%  PMC  400 

0.116 

0.57 

0.535 

Fig.  4 (b) 

Park  [57] 

0.50%  PMC  25 

0.021 

0.63 

0.072 

Fig.  4 (c) 

Al-Fariss  & Pinder  [58] 

waxy  oil  4% 

1.222 

0.77 

3.362 

Fig.  4 (d) 

Al-Fariss  & Pinder  [58] 

waxy  oil  5% 

0.463 

0.87 

3.575 

Fig.  4 (e) 

Chase  & Dachavijit  [59] 

Carbopol  941  1.3% 

0.215 

1.00 

28.46 
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Figure  3:  Comparing  Ellis  (solid)  to  Poiseuille  (dashed)  models  for  converging- 
diverging  rigid  tubes  of  the  given  geometry:  (a)  conic  with  L = 0.15,  Rrn  = 0.01, 
Rm  = 0.02,  pj  = 5000,  po  = 0,  Qe  = 0.187942,  QP  = 0.00452118  (b)  parabolic  with 
L = 0.013,  Rrn  = 0.0017,  Rm  = 0.0025,  p{  = 8000,  pQ  = 6000,  QE  = 3.43037 x 10"5, 
Qp  = 1.16241  x 10~5  (c)  hyperbolic  with  L = 0.03,  Rm  = 0.002,  Rm  = 0.004, 
Pi  = 7000,  po  = 4000,  Qe  = 8.49764  x 10~6,  QP  = 8.0162  x 10'6  (d)  hyperbolic 
cosine  with  L = 0.6,  Rm  = 0.04,  RM  = 0.1,  pt  = 4000,  pQ  = 2000,  Qe  = 1.8855, 
QP  = 0.00184536  (e)  sinusoidal  with  L = 0.55,  Rm  = 0.03,  Rm  = 0.07,  pi  = 15000, 
p0  = 14000,  Qe  = 2.14775,  Qp  = 0.00757797.  All  dimensional  quantities  are 
in  standard  SI  units.  In  all  Eve  sub-figures,  the  vertical  axis  represents  the  axial 
pressure  in  pascals  while  the  horizontal  axis  represents  the  tube  axial  coordinate 
in  meters.  The  rheological  properties  of\he  fluids  are  given  in  Table  2. 


(a)  Conic 


(b)  Parabolic 


Figure  4:  Comparing  Herschel-Bulkley  (solid)  to  Poiscuille  (dashed)  models  for 
converging-diverging  rigid  tubes  of  the  given  geometry:  (a)  conic  with  L = 0.011, 
Rm  = 0.001,  Rm  = 0.0027,  pt  = 1500,  pQ  = 0,  QH  = 4.4286  x 10"4,  QP  = 
2.49962  x 10-6  (b)  parabolic  with  L = 0.65,  Rm  = 0.04,  Rm  = 0.15,  pt  = 7000, 
p0  = 6000,  Qh  = 23.9883,  Qp  = 0.252061  (c)  hyperbolic  with  L = 0.35,  R,rn  = 0.03, 
Rm  = 0.08,  pi  = 10000,  p0  = 5000,  QH  = 0.0625224,  QP  = 0.012097  (d)  hyperbolic 
cosine  with  L = 0.025,  Rm  = 0.0025,  Rm  = 0.005,  Pi  = 8000,  pD  = 5000,  Qp  = 
1.90137  x 10~5,  QP  = 8.13167  x 10"6  (e)  sinusoidal  with  L = 0.05,  Rm  = 0.004, 
Rm  = 0.01,  pi  = 9000,  pa  = 3000,  QH  = 1.84272  x 10"4,  QP  = 2.04623  x 10~4.  All 
dimensional  quantities  are  in  standard  SI  units.  In  all  five  sub-figures,  the  vertical 
axis  represents  the  axial  pressure  in  pascals  while  the  horizontal  axis  represents  the 
tube  axial  coordinate  in  meters.  The  rhjglogical  properties  of  the  fluids  are  given 
in  Table  3. 


Figure  5:  Flow  rate  versus  pressure  drop  for  the  flow  of  a 6.0%  Elvanol  72-51 
solution  modeled  by  an  Ellis  fluid  with  /i0  = 0.185  Pa.s,  a = 2.4  and  /2  = 1025  Pa 
flowing  in  a converging- diverging  conic  tube  with  L = 0.1  m,  Rm  = 0.005  m and 
Rm  = 0.02  m.  The  bulk  rheology  of  the  fluid  is  taken  from  Sadowski  dissertation 

M- 


Figure  6:  Flow  rate  versus  pressure  drop  for  the  flow  of  a guar  gum  solution  of 
0.72%  concentration  modeled  by  an  Ellis  fluid  with  /iG  = 2.672  Pa.s,  a = 3.46  and 
7"i/2  = 9.01  Pa  flowing  in  a converging-diverging  parabolic  tube  with  L = 0.06  m, 
Rm  = 0.005  m and  Rm  = 0.013  m.  The  bulk  rheology  of  the  fluid  is  taken  from 
Balhoff  dissertation  [20]. 
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0.025 


Figure  7:  Flow  rate  versus  pressure  drop  for  the  flow  of  an  aqueous  solution  of 
polyacrylamide  with  0.25%  concentration  modeled  by  an  Ellis  fluid  with  /i0  = 
1.87862  Pa.s,  a = 2.4367  and  ri/2  = 0.5310  Pa  flowing  in  a converging-diverging 
hyperbolic  tube  with  L = 0.025  m,  Rrn  = 0.002  m and  Rm  = 0.006  m.  The  bulk 
rheology  of  the  fluid  is  taken  from  Park  dissertation  [57]. 


Figure  8:  Flow  rate  versus  pressure  drop  for  the  flow  of  a hypothetical  shear  thicken- 
ing Herschel-Bulkley  fluid  with  C = 0.075  Pa.sn,  n = 1.25  and  rD  = 20.0  Pa  flowing 
in  a converging-diverging  hyperbolic  cosine  tube  with  L = 0.75  m,  Rm  = 0.05  m 
and  Rm  = 0.15  m.  The  threshold  yield  pressure  is  about  417  Pa.  Unlike  the 
plots  of  shear-thinning  and  Bingham  fluids,  the  curve  concave  downward  due  to 
the  shear-thickening  nature. 


14 


Figure  9:  Flow  rate  versus  pressure  drop  for  the  flow  of  a hypothetical  shear- 
thinning Herschel-Bulkley  fluid  with  C = 0.673  Pa.sn,  n = 0.54  and  rQ  = 20.0  Pa 
flowing  in  a converging-diverging  sinusoidal  tube  with  L = 0.5  m,  Rm  = 0.015  m 
and  Rm  = 0.06  m.  The  threshold  yield  pressure  is  about  664  Pa. 


(a)  Conic  (b)  Hyperbolic 


(c)  Sinusoidal 

Figure  10:  A graphic  demonstration  of  a sample  of  diverging-converging  tube  ge- 
ometries. 
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Pressure  Drop  (Pa) 


(a) 


(b) 


Figure  11:  Simulated  flow  of  0.6%  Natrosol  - 250H  solution  modeled  by  an  Ellis  fluid 
with  /j,0  = 0.4  Pa.s,  a = 2.168  and  ti/2  = 5.2  Pa  flowing  in  a diverging-converging 
conic  tube,  similar  to  the  one  in  Figure  10  (a),  with  L = 0.06  m,  Rm  = 0.005  m and 
Rm  = 0.015  m.  The  bulk  rheology  of  the  fluid  is  taken  from  Sadowski  dissertation 

N. 


Figure  12:  Simulated  flow  of  a hypothetical  shear-thickening  power-law  fluid  with 
C = 0.75  Pa.sn  and  n = 1.5  flowing  in  a diverging-converging  hyperbolic  tube, 
similar  to  the  one  in  Figure  10  (b),  with  L = 0.07  m,  R,rn  = 0.005  m and  Rm  = 
0.013  m.  The  flow  rate  plot  concave  downward  due  to  the  shear-thickening  rheology. 
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Figure  13:  Simulated  flow  of  an  aqueous  solution  of  Carbopol  941  with  1.0%  con- 
centration modeled  by  a Bingham  fluid  with  C = 0.128  Pa.s  and  rG  = 17.33  Pa 
flowing  in  a diverging-converging  sinusoidal  tube,  similar  to  the  one  in  Figure  10 
(c),  with  L = 0.1  m,  Rm  = 0.009  m and  Rm  = 0.02  m.  The  threshold  yield  pressure 
is  about  260  Pa.  The  bulk  rheology  of  the  fluid  is  taken  from  Chase  and  Dachavijit 
[59]. 
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5 Tests  and  Validations 


Several  consistency  and  validation  tests  have  been  carried  out  to  verify  the  residual- 
based  lubrication  method  and  its  code.  In  the  following  we  outline  these  tests 

• The  sensibility  of  the  method  and  the  integrity  of  the  code  have  been  verified 
by  smooth  convergence  to  a final  solution  with  mesh  refinement  by  employing 
finer  discretization. 

• Tests  have  been  carried  out  on  the  flow  of  non-Newtonian  fluids  through 
straight  cylindrical  tubes  with  constant  radii  as  a special  case  for  the  converging- 
diverging  tubes.  In  all  cases  the  numerical  solutions  converged  to  the  correct 
analytical  solutions  as  given  by  Equations  2 and  4. 

• The  method  and  the  code  have  also  been  verified  by  the  convergence  to  the 
correct  analytical  solutions  for  the  flow  of  Newtonian  fluids,  as  a special 
case  for  the  non-Newtonian  fluids,  in  the  five  prototype  converging-diverging 
geometries.  These  geometries  have  analytical  solutions,  given  in  Table  4, 
that  have  been  derived  and  validated  previously  for  the  flow  of  Newtonian 
fluids  [37,  40].  The  Newtonian  numerical  solutions  have  been  obtained  both 
as  Poiseuille  flow  and  as  Herschel-Bulkley  flow  with  rD  = 0 and  n — 1;  both 
of  these  numerical  solutions  were  identical  to  the  corresponding  analytical 
solution  within  acceptable  numerical  errors. 

• The  convergence  to  the  correct  analytical  solution  has  also  been  verified  for 
the  flow  of  non-Newtonian  power-law  fluids  through  these  five  converging- 
diverging  geometries  using  analytical  expressions  that  have  been  derived  and 
validated  previously  in  [39];  these  analytical  expressions  are  given  in  Table  5. 

A sample  of  these  comparisons  between  the  numerical  and  analytical  solutions 
is  given  in  Figure  14  for  the  flow  of  a typical  power-law  fluid  through  a tube 
with  a conic  geometry.  As  seen,  the  two  solutions  are  virtually  identical 
within  acceptable  numerical  errors. 
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Table  4:  The  equations  describing  the  dependency  of  the  flow  rate  Q on  the  pres- 
sure drop  A p for  the  flow  of  Newtonian  fluids  through  rigid  tubes  with  the  five 
converging- diverging  geometries  of  Table  1.  These  equations  were  derived  and 
validated  previously  in  [40]. 
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Table  5:  The  equations  describing  the  dependency  of  the  flow  rate  Q on  the  pressure 
drop  A p for  the  flow  of  power-law  fluids  in  rigid  tubes  for  the  five  converging- 
diverging  geometries  of  Table  1 where  Fi  is  the  Appcll  hypergeometric  function,  2iq 
is  the  hypergeometric  function  and  Im  is  the  imaginary  part  of  the  given  function. 
These  relations  were  previously  [39]  derived  and  validated. 
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Figure  14:  Comparing  the  numerical  and  analytical  solutions  for  the  flow  of  a 
typical  power-law  fluid  with  C = 0.5  Pa.sn  and  n = 0.75  through  a conically-shaped 
converging- diverging  tube  with  L = 0.15  m,  Rm  = 0.01  m and  Rm  = 0.02  m.  The 
analytical  solution  is  obtained  from  the  first  equation  in  Table  5. 
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6 Computational  Issues  and  Comparison 


Apart  from  yield-stress  fluids,  we  did  not  experience  convergence  difficulties  that 
normally  associate  such  non-linear  solution  schemes  like  those  that  we  confronted 
in  our  previous  investigation  to  the  flow  of  fluids  through  distensible  structures 
[37,  55].  In  almost  all  the  cases  investigated  in  the  current  study,  the  convergence 
was  immediate  as  it  occurred  within  a few  Newton- Raphson  iterations  even  for  the 
high-pressure  flow  regime  cases  with  the  use  of  extreme  boundary  conditions  and 
eccentric  parameters  for  the  flow,  fluids  and  tube  geometry. 

As  for  the  yield-stress  fluids,  there  were  some  convergence  difficulties  and  hence 
to  overcome  these  difficulties  we  used  a number  of  numerical  tricks.  The  most 
effective  of  these  tricks  may  be  to  start  with  solving  the  problem  assuming  zero 
yield-stress.  We  then  use  this  yield-free  solution  as  an  initial  guess  for  solving 
the  initial  problem  with  non-zero  yield-stress.  The  final  flow  solution  in  the  yield- 
stress  fluid  cases  is  conditioned  that  if  it  resulted  in  rD  > rw  on  any  element,  the 
flow  in  the  tube  is  set  to  zero.  It  should  be  remarked  that  for  the  flow  of  yield- 
stress  fluids  through  cylindrical  tubes  with  fixed  radii,  the  yield  condition  where 
the  pressure  drop  just  overcomes  the  yield-stress  and  hence  the  flow  starts  is  given 
by  the  following  relation  [21,  50,  60] 


For  the  high-pressure  flow  regimes,  the  convergence  is  usually  more  difficult 
than  for  the  low-pressure  flow  regimes,  especially  for  yield-stress  fluids.  An  effective 
trick  in  these  cases  is  to  step  up  on  the  pressure  ladder  by  starting  the  process  with 
solving  the  same  problem  but  with  low  pressure  boundary  conditions  where  the 
convergence  is  easy.  This  low-pressure  solution  is  then  used  as  an  initial  guess  for 
the  next  step  with  a higher  boundary  pressure.  On  repeating  this  process  with  the 
use  of  a suitably  divided  pressure  field,  e.g.  100  Pa  increase  per  step,  high-pressure 
flow  regime  problems  can  be  solved  without  convergence  difficulties.  Although  this 
stepping  up  process  incurs  an  extra  computational  cost,  in  most  cases  this  extra 
cost  is  very  low. 

Also  for  the  yield-stress  fluids,  the  access  of  the  flow  regimes  at  the  border 
of  the  threshold  yield  pressure,  which  is  required  for  determining  the  exact  yield 
point,  may  not  be  easy  due  to  the  absence  of  solutions  before  yield  with  possible 
difficulties  in  providing  a sensible  initial  guess.  The  trick  then  is  to  start  from  a 
relatively  high-pressure  flow  regime  where  the  convergence  to  a solution  is  easy  and 


(7) 
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where  the  flow  system  is  expected  to  have  yielded  already.  The  solution  at  this 
regime  is  then  used  as  an  initial  guess  for  stepping  down  on  the  pressure  ladder, 
as  for  stepping  up  process  which  is  outlined  earlier.  By  using  small  pressure  steps, 
e.g.  1 Pa  decrease  per  step,  the  flow  regimes  at  the  very  edge  of  the  threshold 
yield  pressure  can  be  accessed  and  hence  the  exact  threshold  yield  pressure  can  be 
determined. 

Regarding  the  CPU  time  and  memory  requirements,  the  residual-based  method 
has  a very  low  computational  cost.  This  is  partly  due  to  the  nature  of  the  problem 
which  is  related  to  single  tubes  and  hence  it  is  normally  of  limited  size.  The 
discretization  usually  involves  a few  tens  or  at  most  a few  hundreds  of  elements  on 
the  discretization  meshes.  Hence  the  computational  cost  normally  does  not  exceed 
a few  kilobytes  of  memory  and  a few  seconds  of  CPU  time  on  a normal  computer. 

The  advantages  of  the  residual-based  method  over  the  traditional  methods  is 
simplicity,  generality,  ease  of  implementation,  very  good  rate  and  speed  of  con- 
vergence, and  very  low  computational  cost.  Moreover,  the  solutions  obtained  by 
the  residual-based  method  match  in  their  accuracy  any  existing  or  anticipated  an- 
alytical solutions  within  the  given  physical  and  numerical  approximations.  The 
accuracy  of  the  residual-based  method  and  its  convergence  to  the  correct  analyti- 
cal solutions  are  confirmed  in  all  cases  in  which  analytical  solutions  are  available 
such  as  the  limiting  cases  of  straight  geometries  and  Newtonian  fluids,  as  well  as 
the  available  analytical  solutions  for  the  non-Newtonian  power-law  fluids  as  given 
by  the  equations  in  Table  5.  The  main  disadvantage  of  the  residual-based  method 
is  its  one-dimensional  nature  and  hence  it  cannot  be  used  for  obtaining  the  param- 
eters of  the  flow  field  in  dimensions  other  than  the  axial  direction.  The  method 
has  also  other  limitations  related  to  the  physical  assumptions  on  which  the  under- 
lying flow  model  is  based.  However,  most  of  these  limitations  are  shared  by  other 
comparable  methods. 

7 Conclusions 

The  investigation  in  this  paper  led  to  the  proposal,  formulation  and  validation  of  a 
residual-based  lubrication  method  that  can  be  used  to  obtain  the  pressure  field  and 
flow  rate  in  rigid  tubes  with  converging  and  diverging  features,  or  more  generally 
with  a cross  section  that  vary  in  size  and/or  shape  in  the  flow  direction,  for  the 
flow  of  non-Newtonian  fluids. 

The  method  in  its  current  formulation  can  be  applied  to  the  entire  category  of 
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time-independent  non-Newtonian  fluids,  that  is  generalized  Newtonian  fluids.  Two 
time-independent  non-Newtonian  fluids,  namely  Ellis  and  Herschcl-Bulklcy,  have 
been  used  in  this  investigation.  The  method,  in  our  judgement,  has  the  capability 
to  be  extended  to  the  history-dependent  non-Newtonian  fluids,  i.e.  viscoelastic  and 
thixotropic/rheopectic. 

Five  converging- diverging  and  diverging-converging  geometries  have  been  used 
for  test,  validation,  demonstration  and  as  prototypes  for  the  investigation  of  such 
flows.  The  method,  however,  has  a wider  applicability  range  with  regard  to  the 
tube  geometry  as  it  can  be  applied  to  all  geometries  that  vary  in  size  and/or  shape 
in  the  flow  direction  as  long  as  a characteristic  relation,  analytical  or  empirical  or 
numerical,  that  correlates  the  flow  rate  to  the  pressure  drop  over  the  presumably- 
straight  discretized  elements  can  be  found. 

The  residual-based  lubrication  method  was  validated  by  the  convergence  to  the 
correct  analytical  solutions  in  the  special  cases  of  straight  constant-radius  geome- 
tries and  Newtonian  fluids,  as  well  as  the  convergence  to  the  analytical  solutions 
for  the  flow  of  power-law  fluids  through  the  five  converging-diverging  prototype 
geometries.  The  method  has  also  been  endorsed  by  a number  of  qualitatively- 
sensible  tendencies  such  as  the  convergence  to  a final  solution  with  improved  mesh 
refinement. 

The  residual-based  method  has  obvious  advantages  as  compared  to  the  more 
traditional  methods.  These  advantages  include  generality,  ease  of  implementation, 
low  computational  cost,  good  rate  and  speed  of  convergence,  reliability,  and  accu- 
racy. The  main  limitation  of  the  method  is  its  one- dimensional  nature. 
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Nomenclature 

a flow  behavior  index  in  Ellis  model 

7 strain  rate 

/ i generalized  Newtonian  viscosity 

Ho  zero-shear-rate  viscosity 
r shear  stress 

r1/2  shear  stress  when  p = ^ in  Ellis  model 

t0  yield-stress 

tw  shear  stress  at  tube  wall 

C consistency  coefficient  in  Herschcl-Bulklcy  model 
/ flow  continuity  residual  function 

Fi  Appell  hypergeometric  function 
2F\  hypergeometric  function 
J Jacobian  matrix 

L tube  length 

n flow  behavior  index  in  Herschel-Bulkley  model 

N number  of  discretized  tube  nodes 

p pressure 

p pressure  vector 

Pi  inlet  pressure 

p0  outlet  pressure 

A p pressure  drop 

Ap  pressure  perturbation  vector 
Q volumetric  flow  rate 

Qe  numeric  flow  rate  for  Ellis  model 
Qh  numeric  flow  rate  for  Herschel-Bulkley  model 
Qp  numeric  flow  rate  for  Poiseuille  model 

r residual  vector 

R tube  radius 

Rm  minimum  radius  of  converging-diverging  tube 
Rm  maximum  radius  of  converging- diverging  tube 
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x tube  axial  coordinate 
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